Towards a new vision of mesh adaptation methods

and its impact on the simulation of PDEs

Loïc Gouarin et Marc Massot

HPC@Maths team
CMAP / CNRS - Ecole polytechnique

Context

The present work is the result of a team work involving

  • Thomas Bellotti (CNRS, EM2C - Fédé Maths CS)
  • Loïc Gouarin (École polytechnique, CMAP)
  • Josselin Massot (École polytechnique, CMAP)
  • Pierre Matalon (École polytechnique, CMAP)
  • Laurent Séries (École polytechnique, CMAP)
  • Christian Tenaud (CNRS, EM2C - Fédé Maths CS)

Mesh adaptation and error control

Burgers equation - sinus problem

\[ \partial_t u + \partial_x \left ( f(u) \right ) = 0, \quad t \geq 0, \quad x \in \mathbb{R}, \qquad f(u) = \dfrac{u^2}{2}, \]

Consider the Cauchy problem with initial conditions:

\[ u^0(x) = \frac{1}{2} (1+\sin(\pi(x-1))) \quad x \in [-1,1] \]

Adaptive Multiresolution

  • Minimum level \(\underline{\ell}\) and maximum level \(\bar{\ell}\).
  • Cells: \[ C_{\ell, k}:=\prod_{\alpha=1}^d\left[2^{-\ell} k_\alpha, 2^{-\ell}\left(k_\alpha+1\right)\right] \]
  • Finest step: \(\Delta x=2^{-\bar{\ell}}\).
  • Level-wise step: \(\Delta x_{\ell}:=2^{-\ell}=2^{\Delta \ell} \Delta x\).

Wavelets

Decomposition of the solution on a wavelet basis [Daubechies, ’88], [Mallat, ’89] to measure its local regularity. “Practical” approach by [Harten, ’95], [Cohen et al., ’03].

Projection operator

Prediction operator at order \(2 s +1\)

\[ {\hat f}_{\ell+1,2 k}={f}_{\ell, k}+\sum_{\sigma=1}^s \psi_\sigma\left({f}_{\ell, k+\sigma}-{f}_{\ell, k-\sigma}\right) \]

Details are regularity indicator \[ {\mathrm{d}}_{\ell, {k}}:={f}_{\ell, {k}}-{\hat{f}}_{\ell, {k}} \]

Let \(f \in W^{\nu, \infty}\) (neigh. of \(C_{\ell, k}\) ), then \[ \left|{\mathrm{d}}_{\ell, k}\right| \lesssim 2^{-\ell \min (\nu, 2 s +1)}|f|_{W^{\min (\nu, 2 s +1), \infty}} \]

Fast wavelet transform:

means at the finest level can be recast as means at the coarsest level + details \[ \begin{array}{rlr} {f}_{\overline{\ell}} & \Longleftrightarrow & \left({f}_{\underline{\ell}}, {{d}}_{\underline{\ell} +1}, \ldots, {d}_{\bar{\ell}}\right)\\ \end{array} \]

Mesh coarsening (static)

Local regularity of the solution allows to select areas to coarsen

\[ {{f}}_{\bar{\ell}} \rightarrow \left({f}_{\underline{\ell}}, {\mathbf{d}}_{\underline{\ell}+1}, \ldots, {\mathbf{d}}_{\bar{\ell}}\right) \rightarrow \left({f}_{\underline{\ell}}, {\tilde{\mathbf{d}}}_{\underline{\ell}+1}, \ldots, \tilde{{\mathbf{d}}}_{\bar{\ell}}\right) \rightarrow {\tilde{{f}}}_{\bar{\ell}} \] \[ \tilde{{\mathrm{d}}}_{\ell, k}= \begin{cases}0, & \text { if } \left|{\mathbf{d}}_{\ell, k}\right| \leq \epsilon_{\ell}=2^{-d \Delta \ell} \epsilon, \quad \rightarrow \quad\left\|{\mathbf{f}}_{\bar{\ell}}-\tilde{{\mathbf{f}}}_{\bar{\ell}}\right\|_{\ell^p} \lesssim \epsilon \\ {\mathrm{d}}_{\ell, k}, & \text { otherwise} \end{cases} \]

Set a small (below \(\epsilon_{\ell}\)) detail to zero \(\equiv\) erase the cell \(C_{\ell, k}\) from the structure

Examples

Equation

\[ f(x) = 1 - \sqrt{\left| sin \left( \frac{\pi}{2} x \right) \right|} \; \text{for} \; x\in[-1, 1] \]

min level 1
max level 12
ε 10-3
compression rate 96.29%
error 0.00053

Time evolution of PDEs

  • Finite volumes with global time step \(\Delta t = \Lambda(\Delta x)\)
  • Use dynamic mesh refinement

Mesh updated using “old” information at time \(t\) to accommodate the one at time \(t + \Delta t\)

  • Propagation of information : add security cells
  • Formation of singularities : (regularity index: \(\nu =0\), \(\mu = \min(\nu,2 s +1)\)) refine if \[ \left|{\mathbf{d}}_{\ell, k}\right| \geq \epsilon_{\ell}\,2^{d+\mu} \]

Finite volumes / conservation / order

Flux evaluation at interfaces between levels

Using the prediction operator allows to evaluate fluxes at the same level

Finite volume method

We use a Godunov flux for the small hat problem

Numerical Analysis and Modified Equations

Linear scalar transport equation

In this work, we are concerned with the numerical solution of the Cauchy problem associated with the linear scalar conservation law

\[ \partial_t u(t, x)+V \partial_x u(t, x)=0, \quad(t, x) \in \mathbb{R}^{+} \times \mathbb{R} \]

where \(V\) is the transport velocity, taken \(V>0\) without loss of generality. we consider 1 d problems. The extension to \(2\mathrm{~d}\) / \(3\mathrm{~d}\) problems is straightforward and usually done by tensorization [Bellotti2022] and yields analogous conclusions.

The discrete volumes are \[ C_{\ell, k}:=\left[2^{-\ell} k, 2^{-\ell}(k+1)\right], \quad k \in \{ 0,2^{\ell}-1 \}, \] for any \(\ell \in \{ \underline{\ell}, \bar{\ell} \}\). The measure of each cell at level \(\ell\) is \(\Delta x_{\ell}:=2^{-\ell}\) and we shall indicate \(\Delta x:=\Delta x_{\bar{\ell}}\). The cell centers are \(x_{\ell, k}:=\) \(2^{-\ell}(k+1 / 2)\). Finally, we shall indicate \(\Delta \ell:=\bar{\ell}-\ell\), hence \(\Delta x_{\ell}=2^{\Delta \ell} \Delta x\).

Finite Volume scheme

Finite Volume scheme at the finest level of resolution \(\bar{\ell}\) for any cell of indices \(\bar{k} \in \{ 0,2^{\bar{\ell}}-1 \}\). Explicit schemes read:

\[ \mathrm{v}_{\bar{\ell}, \bar{k}}^{n+1}=\mathrm{v}_{\bar{\ell}, \bar{k}}^n-\frac{\Delta t}{\Delta x}\left(\Phi\left(\mathrm{v}_{\bar{\ell}, \bar{k}+1 / 2}^n\right)-\Phi\left(\mathrm{v}_{\bar{\ell}, \bar{k}-1 / 2}^n\right)\right) \]

where we utilize the same linear numerical flux for the left and the right flux (conservativity) \[ \Phi\left(\mathbf{v}_{\bar{\ell}, \bar{k}-1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha \mathbf{v}_{\bar{\ell}, \bar{k}+\alpha}, \quad \Phi\left(\mathbf{v}_{\bar{\ell}, \bar{k}+1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha v_{\bar{\ell}, \bar{k}+1+\alpha} \]

Modified equations

[Carpentier et al 97] or Cauchy-Kowalewski procedure [Harten et al 87]

\[ \partial_t u\left(t^n, x_{\bar{\ell}, \bar{k}}\right)+V \partial_x u\left(t^n, x_{\bar{\ell}, \bar{k}}\right)=\sum_{h=2}^{+\infty} \Delta x^{h-1} \sigma_h \partial_x^h u\left(t^n, x_{\bar{\ell}, \bar{k}}\right) \]

  • Upwind scheme \[ \partial_t u+V \partial_x u=\frac{\Delta x V}{2}(1-\lambda V) \partial_{x x} u+O\left(\Delta x^2\right) \]
  • Lax-Wendroff scheme \[ \partial_t u+V \partial_x u=-\frac{\Delta x^2 V}{6}\left(1-\lambda^2 V^2\right) \partial_x^3 u+O\left(\Delta x^3\right) \]
  • OSMP-3 scheme \[ \partial_t u+V \partial_x u=\frac{\Delta x^3 V}{24}\left(-\lambda^3 V^2+2 \lambda^2 V^2+\lambda V-2\right) \partial_x^4 u+O\left(\Delta x^4\right) \]

How to include MRA I

We introduce the reconstruction operator \(\hat{s}\) instead of \(s\) on the cells \(\left(\bar{\ell}, 2^{\Delta \ell} k+\delta\right)\) for any \(\delta \in \mathbb{Z}\) at the finest level

  • \(\hat{s}=s\) : exact local flux reconstruction [Cohen et al. 2003].
  • \(\hat{s}=0\) but \(s>0\), direct evaluation or naive evaluation [Hovhannisyan et al 2010].

\[ \mathbf{w}_{\bar{\ell}, \bar{k}}^{n+1}=\mathbf{w}_{\bar{\ell}, \bar{k}}^n-\frac{\Delta t}{\Delta x}\left(\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}+1 / 2}^n\right)-\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}-1 / 2}^n\right)\right) \] \[ \Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}-1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha \hat{\hat{\mathbf{w}}}_{\bar{\ell}, \bar{k}+\alpha} \]

How to include MRA II

Let now \((\ell, k) \in S\left(\tilde{\Lambda}^{n+1}\right)\), taking the projection yields the multiresolution scheme

\[ \mathbf{w}_{\ell, k}^{n+1}=\mathbf{w}_{\ell, k}^n-\frac{\Delta t}{\Delta x_{\ell}}\left(\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell}(k+1)+1 / 2}^n\right)-\Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell} k-1 / 2}^n\right)\right) \]

\[ \Phi\left(\hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell} k-1 / 2}\right):=V \sum_{\alpha=\underline{\alpha}}^{\bar{\alpha}} \phi_\alpha \hat{\hat{\mathbf{w}}}_{\bar{\ell}, 2^{\Delta \ell} k+\alpha} \]

Some information is loss because of the averaging procedure: two different schemes we can consider for the computation of the modified equations.

Theorem

The local truncation error of the reference Finite Volume scheme and the one of the adaptive Finite Volume scheme are the same up to order \(2\hat{s}+1\) included.

Modified equations including MRA

This result establishes at which order the modified equations of the reference scheme are perturbed by the introduction of the adaptive scheme. However, it does not characterize the terms in the modified equations above order \(2\hat{s}+1\) in \(\Delta x\) (symbolic computations).

  • Upwind scheme (original) \[ {\color{blue}{% \partial_t u+V \partial_x u=\frac{\Delta x V}{2}(1-\lambda V) \partial_{x x} u+O\left(\Delta x^2\right) }} \]
  • Upwind scheme with MRA \[ \begin{array}{lr} \partial_t u+V \partial_x u=\frac{\Delta x V}{2}\left(2^{\Delta \ell}-\lambda V\right) \partial_{x x} u+O\left(\Delta x^2\right), & \text { for } \hat{s}=0 \\ \partial_t u+V \partial_x u=\frac{\Delta x V}{2}(1-\lambda V) \partial_{x x} u-\frac{\Delta x^2 V}{6}\left(1-\lambda^2 V^2\right) \partial_x^3 u+& \\ \ \ \ \ \ \ \ \ \ \ \frac{\Delta x^3 V}{24}\left(-3 \Delta \ell 2^{2 \Delta \ell}+2^{2 \Delta \ell}-\lambda^3 V^3\right) \partial_x^4 u+O\left(\Delta x^4\right), & \text { for } \hat{s}=1 \end{array} \]

Modified equations including MRA

This result establishes at which order the modified equations of the reference scheme are perturbed by the introduction of the adaptive scheme. However, it does not characterize the terms in the modified equations above order \(2\hat{s}+1\) in \(\Delta x\) (symbolic computations).

  • Lax-Wendroff scheme (original) \[ {\color{blue}{\partial_t u+V \partial_x u=-\frac{\Delta x^2 V}{6}\left(1-\lambda^2 V^2\right) \partial_x^3 u+O\left(\Delta x^3\right)}} \]
  • Lax-Wendroff scheme with MRA \[ \begin{array}{lr} \partial_t u+V \partial_x u=\frac{\Delta x \lambda V^2}{2}\left(2^{\Delta \ell}-1\right) \partial_{x x} u+O\left(\Delta x^2\right), & \text { for } \hat{s}=0 \\ \partial_t u+V \partial_x u=-\frac{\Delta x^2 V}{6}\left(1-\lambda^2 V^2\right) \partial_x^3 u+&\\ \ \ \ \ \ \ \ \ \ \ \frac{\Delta x^3 \lambda V^2}{24}\left(-3 \Delta \ell 2^{2 \Delta \ell}+2^{2 \Delta \ell}-\lambda^2 V^2\right) \partial_x^4 u+O\left(\Delta x^4\right), & \text { for } \hat{s}=1 \end{array} \]

Theoretical results on the global error

Theorem 2

Assume that

  • The reference scheme satisfies the restricted stability condition \(\|E\| \leq 1\)
  • The Harten-like scheme satisfies the restricted stability condition \(\left\|\bar{E}_{\Lambda}\right\| \leq 1\) for any \(\Lambda\).

Then, for smooth solution, in the limit \(\Delta x \rightarrow 0\) (i.e. \(\bar{\ell} \rightarrow+\infty\) ) and for \(\Delta \underline{\ell}=\bar{\ell}-\underline{\ell}\) kept fixed, we have the error estimate

\[ \left\|\mathbf{v}_{\bar{\ell}}^n-\mathbf{w}_{\bar{\ell}}^n\right\| \leq C_{t r} t^n \Delta x^{2 \hat{s}+1}+C_{m r} \frac{t^n}{\lambda \Delta x} \epsilon \]

where \(C_{t r}=C_{t r}\left(\bar{\ell}-\underline{\ell},\left(\phi_\alpha\right)_\alpha, \lambda, \hat{s}, V\right)\) and \(C_{m r}=C_{m r}\left(\bar{\ell}-\underline{\ell},\left(\phi_\alpha\right)_\alpha, \lambda, \hat{s}, s, V\right)\). \[ \left\|\mathbf{u}_{\bar{\ell}}^n-\mathbf{w}_{\bar{\ell}}^n\right\| \leq C_{r e f} t^n \Delta x^\theta+C_{t r} t^n \Delta x^{2 \hat{s}+1}+C_{m r} \frac{t^n}{\lambda \Delta x} \epsilon \]

Comments on theorem

  • The error estimate contains three contributions: the discretization error of the reference scheme, the perturbation error between the reference and the adaptive scheme, and the thresholding error coming from the multiresolution
  • The constant \(C_{\mathrm{tr}}\) generally grows exponentially with \(\bar{\ell}-\underline{\ell}\), sometimes also involving linear terms, i.e. \(\hat{s}=1\). We have the following cases:
    • \(\theta<2 \hat{s}+1\). The error of the reference scheme dominates the perturbation introduced by the adaptive scheme \(\left\|\mathbf{u}_{\bar{\ell}}^N-\mathbf{w}_{\bar{\ell}}^N\right\| \leq C_{\mathrm{ref}} T \Delta x^\theta+C_{\mathrm{mr}} \frac{T}{\lambda \Delta x} \epsilon\). A thresholding error of the same order as the reference error \(\epsilon \sim \Delta x^{\theta+1}\).
    • \(\theta=2 \hat{s}+1\). The error of the reference scheme and the perturbation order are comparable (first example!) We have \(\left\|\mathbf{u}_{\bar{\ell}}^N-\mathbf{w}_{\bar{\ell}}^N\right\| \leq\left(C_{\mathrm{ref}}+C_{\mathrm{tr}}\right) T \Delta x^\theta+C_{\mathrm{mr}} \frac{T}{\lambda \Delta x} \epsilon\).
    • \(\theta>2 \hat{s}+1\). The perturbation introduced by the adaptive scheme dominates the error of the reference scheme. Therefore, multiresolution introduces a large perturbation that yields a different convergence rate. We have \(\left\|\mathbf{u}_{\bar{\ell}}^N-\mathbf{w}_{\bar{\ell}}^N\right\| \leq C_{\mathrm{tr}} T \Delta x^{2 \hat{s}+1}+C_{\mathrm{mr}} \frac{T}{\lambda \Delta x} \epsilon\), thus \(\epsilon \sim \Delta x^{2 \hat{s}+2}\) (AMR !)

How to compute fluxes at the finest level

How to compute fluxes at the finest level

How to compute fluxes at the finest level

How to compute fluxes at the finest level

Burgers results

Burgers results (Error for scheme order 1)






Burgers results (Error for scheme order 1)






Burgers results (Error for scheme order 1)






Burgers results (MR solution and level for scheme order 1)

Burgers results (MR solution and error for scheme order 1)

Burgers results (MR+MLF sol. and level for scheme order 1)

Burgers results (MR+MLF sol. and error for scheme order 1)

Burgers results (Error for scheme order 1)






Burgers results (MR solution and level for scheme order 1)

Burgers results (MR solution and error for scheme order 1)

Burgers results (MR+MLF sol. and level for scheme order 1)

Burgers results (MR+MLF sol. and error for scheme order 1)

Burgers results (Error for scheme order 3)






Burgers results (Error for scheme order 3)






Burgers results (MR solution and level for scheme order 3)

Burgers results (MR solution and error for scheme order 3)

Burgers results (MR+MLF sol. and level for scheme order 3)

Burgers results (MR+MLF sol. and error for scheme order 3)

Burgers results (MR VS MR+MLF error for scheme order 3)

Burgers 2D results (MR+MLF solution order 1)

Euler 2D results (MR+MLF solution order 3 (OSMP scheme))

Euler 2D results (MR+MLF solution order 3 (OSMP scheme))

Conclusion

  • Analysis of mesh adaptation completed for convection - extension to reaction and diffusion
    • Link with AMR important in terms of accuracy of the simulations
  • Conducting simulation with strict error control requires adpative multiresolution
    • Data structure is a key issue - see presentation of L. Gouarin
    • New generation of coupled splitting / ImEx schemes in samurai
  • Link through the code and developper team to various applications

samurai

Patched-based representation

Advantages

  • compact data structure
  • rectangular zones work well for tiling and optimizing caches
  • efficient to work with neighbors and the different levels

Disadvantages

  • requires more cells than necessary
    • grid hierarchy: overlapping
    • larger refinement zone

Cell-based representation

Advantages

  • no grids hierarchy
  • keep only needed cells

Disadvantages

  • more difficult to find the neighbors
  • more difficult to work with the different levels
  • tree-like structure
  • more difficult to add ghost cells hierarchy

Adaptive methods

AMR

Advantages

  • based on heuristic criteria (gradient, second derivative, …)
  • easy to implement

Disadvantages

  • physical problem understanding is important
  • add potentially more cells that needed

MRA (see the presentation of M. Massot)

Advantages

  • based on wavelet decomposition
  • error control
  • independent of the physical problem
  • only needed cells

Disadvantages

  • difficult to implement with the current data structures
  • can be costly

Can we have a data structure that takes the advantages of both and well suited to implement various adaptive mesh refinement methods ?

Samurai

Design principles

  • Compress the mesh according to the level-wise spatial connectivity along each Cartesian axis.
  • Achieve fast look-up for a cell into the structure, especially for parents and neighbors.
  • Maximize the memory contiguity of the stored data to enable caching and vectorization.
  • Manage and facilitate complex operations between meshes (numerical schemes, MRA, …).

An overview of the data structure

[ start , end [ @ offset

An overview of the data structure

  • Level 0
    • x: [0, 4[, [0, 1[, [3, 4[, [0, 1[, [3, 4[, [0, 3[
    • y: [0, 4[
    • y-offset: [0, 1, 3, 5, 6]
  • Level 1
    • x: [2, 6[, [2, 6[, [2, 4[, [5, 6[, [2, 6[, [6, 8[, [6, 7[
    • y: [2, 8[
    • y-offset = [0, 1, 2, 4, 5, 6, 7]
  • Level 2
    • x: [8, 10[, [8, 10[, [14, 16[, [14, 16[
    • y: [8, 10[, [14, 16[
    • y-offset: [0, 1, 2, 3, 4]

Identify the different types of cells

Identify the different types of cells

Identify the different types of cells

Identify the different types of cells

Algebra of sets

\

=

Algebra of sets

The search of an admissible set is recursive. The algorithm starts from the last dimension (y in 2d, z in 3d,…).

The available operators in samurai are for now

  • the intersection of sets,
  • the union of sets,
  • the difference between two sets,
  • the translation of a set,
  • the projection of a set.

Usage example: find jumps

auto jump_set = intersection(
                    translate(mesh[1], {1}).on(0),
                    mesh[0]
                );

Usage example: the projection operator

auto proj_set = intersection(mesh[level], mesh[level + 1])
          .on(level);

proj_set([&](auto i)
{
    u(level, i) = 0.5*(u(level+1, 2*i) + u(level+1, 2*i+1));
});

Compression rates

Compression rates

Level Num. of cells p4est samurai (leaves) samurai (all) ratio
\(9\) 66379 2.57 Mb 33.68 Kb 121 Kb 21.24
\(10\) 263767 10.25 Mb 66.64 Kb 236.8 Kb 43.28
\(11\) 1051747 40.96 Mb 132.36 Kb 467.24 Kb 87.66
\(12\) 4200559 163.75 Mb 263.6 Kb 927 Kb 176.64
\(13\) 16789627 654.86 Mb 525.9 Kb 1.85 Mb 353.98
\(14\) 67133575 2.61 Gb 1.05 Mb 3.68 Mb 709.24

Other features

  • Loop algorithms over the levels and the cells
  • Simplified access operator
  • Reconstruct a value anywhere even if the cell doesn’t exist
  • Helper classes to construct complex meshes
  • Helper classes to construct schemes for explicit and implicit usage
  • Helper classes to construct N-D operators and expressions using xtensor or Eigen
  • MPI support
  • HDF5 support

Roadmap

People involved

  • 4 core developers
  • users and developers communities
  • Various collaborations and research programs

Real-World Applications

  • Interfacial flow simulation
  • Simulation analysis on the Hydrogen risk
  • Plasma dynamics
  • Battery simulation
  • Combustion modeling